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HIV dynamic studies have contributed significantly to the under- 
standing of HIV pathogenesis and antiviral treatment strategies for 
AIDS patients. Establishing the relationship of virologic responses 
with clinical factors and covariates during long-term antiretroviral 
(ARV) therapy is important to the development of effective treat- 
ments. Medication adherence is an important predictor of the effec- 
tiveness of ARV treatment, but an appropriate determinant of ad- 
herence rate based on medication event monitoring system (MEMS) 
data is critical to predict virologic outcomes. The primary objective 
of this paper is to investigate the effects of a number of summary 
determinants of MEMS adherence rates on virologic response mea- 
sured repeatedly over time in HIV-infected patients. We developed 
a mechanism-based differential equation model with consideration 
of drug adherence, interacted by virus susceptibility to drug and 
baseline characteristics, to characterize the long-term virologic re- 
sponses after initiation of therapy. This model fully integrates viral 
load, MEMS adherence, drug resistance and baseline covariates into 
the data analysis. In this study we employed the proposed model and 
associated Bayesian nonlinear mixed-effects modeling approach to as- 
sess how to efficiently use the MEMS adherence data for prediction of 
virologic response, and to evaluate the predicting power of each sum- 
mary metric of the MEMS adherence rates. In particular, we intend 
to address the questions: (i) how to summarize the MEMS adherence 
data for efficient prediction of virologic response after accounting for 
potential confounding factors such as drug resistance and covariates, 



Received May 2009; revised September 2009. 

^Supported in part by NIAID/NIH Grant AI080338 and MSP/NSA Grant H98230-09- 
1-0053 to Y. Huang, and NIH Grants AI50020, AI078498, AI078842 and AI087135 to H. 
Wu. 

Key words and phrases. Bayesian mixed-effects models, confounding factors, HIV dy- 
namics, longitudinal data, MEMS adherence assessment, time- varying drug efficacy, virus 
resistance. 

This is an electronic reprint of the original article published by the 
Institute of Mathematical Statistics in The Annals of Applied Statistics, 
2011, Vol. 5, No. 1, 551-577. This reprint differs from the original in pagination 
and typographic detail. 



1 



2 



HUANG ET AL. 



and (ii) how to evaluate treatment effect of baseline characteristics 
interacted with adherence and other clinical factors. The approach 
is applied to an AIDS clinical trial involving 31 patients who had 
available data as required for the proposed model. Results demon- 
strate that the appropriate determinants of MEMS adherence rates 
are important in order to more efficiently predict virologic response, 
and investigations of adherence to ARV treatment would benefit from 
measuring not only adherence rate but also its summary metric as- 
sessment. Our study also shows that the mechanism-based dynamic 
model is powerful and efi'ective to establish a relationship of virologic 
responses with medication adherence, virus resistance to drug and 
baseline covariates. 

1. Introduction. The revolution in human immunodeficiency virus (HIV) 
treatment has brought diagnostic tests that can accurately measure levels of 
HIV in blood. Resulting data show (plasma) viral load (HIV-1 RNA copies or 
RNA copies) to be an important predictor of the risk of progression to AIDS. 
The antiretroviral (ARV) agents, which include potent protease inhibitors 
(Pis) are, however, not a cure for HIV infection. While many patients benefit 
from ARV treatment, others do not benefit or only experience a temporary 
benefit. There are several reasons why treatment fails, of which poor patient 
adherence to ARV therapy is a leading factor [Ickovics and Meisler (1997); 
Paterson et al. (2000)]. Thus, assessment of adherence within AIDS clinical 
trials is a critical component of the successful evaluation of therapy out- 
comes. Maintaining adherence may be particularly difficult when the drug 
regimen is complex or side-effects are severe, as is often the case for current 
HIV therapy [Ickovics and Meisler (1997)]. 

The measurement of adherence remains problematic; a standard definition 
of adherence and completely reliable measures of adherence are lacking. 
Nevertheless, there has been substantial progress in both of these areas in 
the past few years. First, it appears that higher levels of adherence are 
needed for HIV disease than other diseases to achieve the desired therapeu- 
tic benefit [Paterson et al. (2000)]. Second, better appreciation of the value 
and limitations of different adherence measurements has been addressed 
[Bova et al. (2005)]. In AIDS clinical trials adherence to medication regimen 
is currently measured by two methods: by use of questionnaires (patient 
self-reporting and/or face-to- face interview) and by use of electronic com- 
pliance monitoring (Medication Event Monitoring System [MEMS]) caps. 
The MEMS is often used as an objective adherence measure. It consists of 
a computer chip in the cap of a medication bottle that records each time 
the bottle is opened. The results can be downloaded, printed out and ana- 
lyzed. It demonstrates that medication-taking patterns are highly variable 
among patients [Kastrissios et al. (1998)] and that they often give a more 
precise measure of adherence than self-report [Arnsten et al. (2001)]. How- 
ever, MEMS data are also subject to error and are not widely available in 
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the clinical setting. Adherence assessment by self-report is usually evalu- 
ated by a patient's ability to recall their medication dosing during a specific 
time interval. Finally, it is important to note that the measurement of viral 
load levels is of special utility as an indirect measure of adherence in HIV 
therapeutics. It has been argued that this is not an adherence measure be- 
cause other factors may influence viral load (drug resistance, etc.). However, 
there is a tight correlation between viral load and adherence [Paterson et al. 
(2000); Haubrich et al. (1999)]. 

Viral dynamic models can be formulated through ordinary differential 
equations (ODE) , but there has been only limited development of statistical 
methodologies for assessing their agreement with observed data. Currently 
there also are substantial knowledge gaps between theoretical HIV dynam- 
ics and the role of many clinical factors. In developing long-term dynamic 
modeling, this paper will address these problems by utilizing time-specific 
information, such as drug adherence and susceptibility factors, on the bio- 
logical mechanism of HIV dynamics to achieve more realistic and accurate 
characterization of the relationship between clinical/drug factors and viro- 
logic response. Several studies [Arnsten et al. (2001); Levine et al. (2006)] 
investigated the association between virologic responses and adherence as- 
sessed by both questionnaire and MEMS data. The results indicated that 
the MEMS cap adherence data may not be correlated better to virologic re- 
sponse compared to the questionnaire adherence data unless the MEMS cap 
data are summarized in an appropriate way. Further, Huang et al. (2008), 
Labbe and Verotta (2006), Liu et al. (2007) and Vrijens et al. (2005) modeled 
the relationship between virologic response and adherence rate using ques- 
tionnaire data and MEMS data averaged by each interval between study 
visits or weekly basis, but no significant differences were found in predicting 
virologic response. Along with this line, this paper will investigate different 
determinants of the adherence rate based on MEMS data from an AIDS clin- 
ical trial study [Hammer et al. (2002)] and compare their performance for 
predicting a virological response. We employed the proposed mechanism- 
based dynamic model to assess how to efficiently use the adherence data 
based on MEMS to predict virological response. In particular, we intend 
to address the questions (i) how to summarize the MEMS adherence data 
for efficient prediction of virological response after accounting for potential 
confounding factors such as drug resistance and baseline covariates, and (ii) 
how to evaluate treatment effect of baseline characteristics interacted with 
MEMS adherence and other clinical factors. 

The purpose of this paper is to describe a reparameterized ODE dynamic 
model (with identifiable parameters) which fully integrates viral load, medi- 
cation adherence, drug resistance and baseline covariates data from an AIDS 
clinical trial study into the analysis. Thus, our dynamic model will be able 
to characterize sustained suppression or resurgence of the virus as arising 
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from intrinsic viral dynamics, and/or influenced by factors such as drug sus- 
ceptibility and adherence during the treatment period of the clinical trial. 
The Bayesian nonlinear mixed-effects (BNLME) modeling approach [David- 
ian and Giltinan (1995)] is employed to estimate dynamic parameters and 
identify significant clinical factors and/or covariates on virologic response 
to ARV treatment. The rest of this article is organized as follows. Section 
2 introduces reparameterized viral dynamic models with time-varying drug 
efficacy which incorporates the effects of drug adherence, drug resistance 
and baseline covariates, and briefly describes the BNLME modeling ap- 
proach, implemented via Markov chain Monte Carlo (MCMC) procedures, 
followed by defining the deviance information criterion (DIG) for compari- 
son of models. In Section 3 we summarize the motivating data set from an 
AIDS clinical trial study including the data of plasma viral load, medication 
adherence from MEMS cap, drug resistance and baseline covariates; the pro- 
posed methodology is applied to these data and the results are presented. 
The method is evaluated via a simulation study in Section 4. Finally, we 
conclude the article with some discussions in Section 5. 

2. HIV dynamic mechanism-based ODE models and statistical approaches. 

This section aims to introduce long-term viral dynamic models based on 
a system of ODE with time-varying coefficients but without closed-form 
solutions, and to investigate associated methodologies to demonstrate the 
application of these models to an AIDS clinical trial study. Long-term vi- 
ral dynamic models can be used to describe the interaction between cells 
susceptible to target cells (T), infected cells (T*) and free virus (V) by 
considering time- varying drug efficacy [Huang and Wu (2006); Huang et al. 
(2006)]. These three compartments (variables) are described as follows. 

HIV virions (V) will infect target cells (T) and turn them into infected 
cells (T*) at an infection rate k. Due to the intervention of antiviral drugs, 
we assume that drugs reduce the infection rate in the infected cells (T*) by 
1 ~ 7(^) [0 — 7(^) — The infected cells will die at rate 5 after producing 
an average of N virions per cell during their lifetimes, and free virions are 
removed from the system at rate c. In addition to the dynamics describing 
virus infection, we have to specify the dynamics of the uninfected cell pop- 
ulation. The simplest assumption is that uninfected cells are produced at 
a constant rate A at which new T cells are generated from sources within 
the body, such as the thymus and die at a rate dx- Thus, the HIV dynamic 
model, after initiation of antiviral therapy, can be written as 

^Tit) = A - drm - [1 - ^it)]kTit)Vit), 
(1) ^T*{t) = [l-j{t)]kT{t)V{t)-6T*{t), 
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—V(t) = N6T*(t)-cV(t) 
at 



where the time- varying parameter 7(t) (as defined below) quantifies the 
time- varying drug efficacy. If the regimen is not 100% effective [i.e., < 
7(t) < 1], the system of ODE cannot be solved analytically. The solutions 
to (1) then have to be evaluated numerically. When ^(t) = 70 (an unknown 
constant), the model (1) becomes the model developed by Perelson and 
Nelson (1999). In particular, when ^(t) = (the drug has no effect), the 
model (1) reduces to the model in the publications [Bonhoeffer et al. (1997); 
Nowak et al. (1995, 1997, 2000); Stafford et al. (2000)]; while -f{t) = 1 (the 
drug is 100% effective), the model (1) reverts to the model discussed by 
Nowak and May (2000) and Perelson and Nelson (1999). 

However, it is challenging to estimate all the parameters in the model (1) 
and to conduct inference because the model (1) is not a priori identifiable 
(i.e., multiple sets of parameters obtain identical fits to the data), given 
only viral load measurements [Cobelli et al. (1979)]. To obtain a model 
with a priori identifiable parameters [Labbe and Verttoa (2006)], this paper 
investigates mechanism-based reparameterized ODE models to quantify the 
long-term viral dynamics with ARV treatment and the associated statistical 
methods for model fitting. 

2.1. Reparameterized model with time-varying drug efficacy. Following 
the studies [Perelson and Nelson (1999); Nowak and May (2000); Labbe and 
Verttoa (2006)], we reparameterize the model (1) using the rescaled variables 
f{t) = {dT/X)T{t), f*{t) = {5/X)T*{t), V{t) = {k/dT)V{t). These yield the 
rescaled version as follows: 



where R = kN\/{cdT) represents the basic reproductive ratio for the virus, 
defined as the number of newly infected cells that arise from any one infected 
cell when almost all cells are uninfected [Nowak and May (2000)]. Note that 
the rescaled model (2) has fewer parameters than the 'original' model (1). 
The identifiability of the model (2) is guaranteed [Cobelli et al. (1979); 
Labbe and Verttoa (2006)] and parameters of the model can be uniquely 
identified. If i? < 1, then the virus will not spread, since every infected cell 
will on average produce less than one infected cell. If, on the other hand, 
i? > 1, then every infected cell will on average produce more than one newly 



(2) 




-jjT = dT{i - m - [1 - 7(t)]T(t)y(t)} 

jjT*{t) = 6{[i - mmvit) - f*{t)}, 

±±V{t) = c{Rf*{t)-V{t)}, 
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infected cell and the virus will proliferate. For the HIV virus to persist in 
the host, infected cells must produce at least one secondary infection, and 
R must be greater than unity [Nowak and May (2000)]. 

Assuming steady state before the beginning of drug therapy, initial con- 
ditions for the model can riow be expressed as simple functions^ of the initial 
conditions for viral load (Vo): Tq = 1/(1 + Vo),T^ = Vb/(1 + Vq), Vb = ^(0) 
[Cobelli et al. (1979); Labbe and Verttoa (2006)]. The assumption of initial 
steady state is necessary to guarantee identifiable (none of the models re- 
ported or referenced here is identifiable if the initial states are unknown), 
and is often justified by the clinical trial protocol. For example, in ACTG 
398, individual patients were taken off the drug before the initiation of the 
new therapy (washout period to eliminate the effect of previously adminis- 
tered drugs and to guarantee that all individuals started from steady-state 
conditions). Finally, viral load [V^(t)] in model (1) is related to an equation 
output of viral load amount [V{t)] in model (2) as follows: V{t) = pV{t), 
where p, which is equivalent to a volume of distribution of pharmacokinet- 
ics, is a viral load scaling (proportionality) factor (10,000 copies/ml) to be 
estimated from the data [Nowak and May (2000)]. The set of ODE (2) wiU 
be used to construct the BNLME model. 

2.2. Time-varying drug efficacy model. Within the population of HIV 
virions in a human host, there is likely to be genetic diversity and corre- 
sponding diversity in susceptibility to the various ARV agents. In clinical 
practice, genotypic or phenotypic tests can be performed to determine the 
sensitivity of HIV-1 RNA to ARV agents before a treatment regimen is 
selected. Here we use the phenotypic marker, /C50 [Molla et al. (1996)], to 
quantify agent-specific drug susceptibility. Because experimental data track- 
ing development of resistance suggest that the resistant fraction of the viral 
population grows exponentially, we propose a model of log-linear function 
to approximate the within- host changes over time in /C50 as follows: 



where and Sj are respective values of IC^oit) at baseline and time point 
at which the resistant mutations dominate. In our study, t/ is the time 
of virologic failure which is observed from clinical studies. If Sf = Sq, no 
new drug resistant mutation is developed during treatment. Although more 
complicated models for median inhibitory concentration have been proposed 
based on the frequencies of resistant mutations and cross-resistance patterns 
[Wainberg et al. (1996); Bonhoeffer, Lipsitch and Levin (1997)], in clinical 
studies or clinical practice it is common to collect /C50 values only at base- 
line and failure time as designed in ACTG 5055 [Acosta et al. (2004)] and 




for <t <tf 
for t>tf, 
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ACTG 398 [Hammer et al. (2002); Pfister et al. (2003)]. Thus, given that 
/C50 is only measured at basehne and at the time of treatment failure, this 
function may serve as a good approximation in terms of data availability. 

Poor adherence to a treatment regimen is one of the major causes of 
treatment failure [Ickovics and Meisler (1997)]. The following model is used 
to represent adherence for a time interval Tk <t < T^+i- 



Ait) 



1, if all doses are taken in (T^,, T^+i], 
rk, if 100rfc% doses are taken in {Tk, Tk+i], 

where < < 1, with indicating the adherence rate computed for each 
assessment interval {Tk,Tk+i] between study visits based on the question- 
naire or MEMS data; Tj. denotes the kth adherence assessment time. 

In most viral dynamic studies, investigators assumed that either drug effi- 
cacy was constant over treatment time [Perelson and Nelson (1999); Wu and 
Ding (1999)] or antiviral regimens had perfect effect in blocking viral repli- 
cation [Ho et al. (1995); Perelson et al. (1996)]. However, the drug efficacy 
may change as concentrations of ARV drugs and other factors (e.g., drug 
resistance) vary during treatment. A simple pharmacodynamic sigmoidal 
-E'max model for dose-effect relationship is [Gabrielsson and Weiner (2000)] 

(3) E = 

EC^o + C 

where -Emax is the maximal effect that can be achieved, C is the drug concen- 
tration, and EC50 is the drug concentration that induced an effect equivalent 
to 50% of the maximal effect. Many different variations of the -Emax model 
have been developed by pharmacologists to model pharmacodynamic effects. 
Ema.x models include the sigmoid -Emax model, the ordinary -Emax model and 
composite E'max models [Gabrielsson and Weiner (2000); Davidian and Gilti- 
nan (1995)]. The ordinary E^max model describes agonistic and antagonistic 
(inhibitory) effects of a drug, the sigmoid E'max model is more flexible for 
the steepness or curvature of the response-concentration curve compared to 
the ordinary Emax model, and composite E^max models are used for multi- 
ple drug effects. More detailed discussions on Smax models can be found in 
the book by Gabrielsson and Weiner (2000) and the paper by Huang et al. 
(2003). Here we employ the following modified E'max model to represent the 
time- varying drug efficacy for two ARV agents within a class, 

A^{t)/ICUt)+A2{t)/ICl,{t) 



cP + Ai{t)/ICUt) + A2{t)/IClo{ty 



where Ak{t) and /C5o(t) {k = 1,2) are the adherence profile of the drug as 
measured by MEMS data and the time-course of median inhibitory con- 
centrations for the two drugs, respectively; (f) = exp(/3o + f^iwi + /32W2)] wi 
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and W2 are observed baseline viral load and CD4 cell count, respectively; 
f3 = {(3o, Pi, ^2)'^ are unknown covariate effect parameters to be estimated 
from clinical data. If /3i = /32 = (without considering effect of covariates), 
<j) = exp(/3o) can be used to quantify the conversion between in vitro and in 
vivo /C50 which is the case discussed by Huang et al. (2003). If j{t) = 1, 
the drug is 100% effective, whereas if 7(t) = 0, the drug has no effect. Note 
that, if Ak{t), /C5o(t), wi and W2 are measured or obtained from a clinical 
study and (3 can be estimated from clinical data, then the time-varying drug 
efficacy 'y{t) can be estimated for the whole period of ARV treatment. 

2.3. Bayesian modeling approaches. A number of studies investigated 
various statistical methods, including Bayesian approaches, to fit viral dy- 
namic models and to predict virological responses [Han et al. (2002); Huang 
et al. (2006); Perelson et al. (1996); Wu et al. (1998); Wu and Ding (1999)]. 
The Bayesian approach to viral dynamic modeling is particularly appealing 
from a biological perspective, as it allows informative prior distributions to 
be incorporated. From a statistical estimation point of view, a Bayesian ap- 
proach is preferable because of the difficulties which are often encountered 
from a classical approach when models involve the large numbers of param- 
eters, and complex nonlinearity of the subject-specific models. A Bayesian 
nonlinear mixed-effects (BNLME) model allows us to incorporate prior infor- 
mation at the population level into the estimates of dynamic parameters for 
individual subjects. We briefly summarize the main concepts in the Bayesian 
approach to inference and the presentation is, of course, far from exhaus- 
tive [Davidian and Giltinan (1995); Gelfand and Smith (1990); Huang et al. 
(2006); Wakefield et al. (1994); Wakefield (1996)]. 

In reference to the model (2), we denote the number of subjects by 
n and the number of measurements on the ith subject by mj. Let ^l = 
(log c, log 5, log (ir, logp, log -R, ^0, /3i, /32)^, = {^i, « = 1, • • • , n}, 0i = (log q, 
log (5i, log c^Ti, log Pi, log log (^i)'^ and Y = {yij,i = 1, . . . ,n; j = 1, . . . ,mi}. 
Let fij{Oi,tj) = logiQ{V{6i,tj)), where V{6i,tj) is proportional to the nu- 
merical solution of V{t) in the differential equations (2) for the ith subject 
at time tj. Let yij{t) and ej(tj) denote the repeated measurements of viral 
load in log^g scale and a measurement error with mean zero, respectively. 
Note that log-transformation of dynamic parameters and viral load is used 
to make sure that estimates of dynamic parameters are positive and to sta- 
bilize the variance and convergence, respectively. The BNLME model can 
be written in the following three levels [Gelfand and Smith (1990); Davidian 
and Giltinan (1995); Huang and Wu (2006); Wakefield (1996)]. 

Level 1. Within-subject variation: 



(5) 
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where = (yniti), . . ■ , yim.i'tmjV , = {fii{OiM), ■ ■ ■ , fim,iOi,tm,))'^ , 

Level 2. Between-subject variation: 

(6) 6>, = W,/i + b„ [bi|5]]~AA(0,I]), 

where bj are random effects with mean zero. It is noteworthy that, for 
(3 = {Po, Pi, ^2)"^ , no log-transformation is required as they are not nec- 
essarily positive. Wj = (Ig, Jij, J2i), where Ig is an identity matrix and 
3 si = (0, 0, 0, 0, 0, Wsi)'^ (s = 1, 2; 2 = 1, 2, . . . , n) are 6x1 vector, with wu 
and W2i being (standardized) individual baseline viral load (in log^Q scale) 
and CD4 cell count, respectively. For (3 = (/3o, /32)"^, we are only interested 
in estimating them at population level. Thus, the individual parameter (pi 
is related to them as follows, log^j = /3o + PiWu + /32W2i + biQ, where biQ is 
the last element of bj (i = 1, 2, . . . , n). 
Level 3. Hyperprior distributions: 

(7) a-^^Ga{a,b), ~ A/'(r/, A), ~ I^i(n, i/), 

where the mutually independent Gamma (Ga), Normal (Af) and Wishart 
( Wi) prior distributions are chosen to facilitate computations [Davidian and 
Giltinan (1995)]. The hyper-parameters a, b, rj, A, and z/ can be determined 
from previous studies and literature. 

The Bayesian approach is developed in the presence of observations whose 
value is initially uncertain and described through a probability distribution, 
which depends on some parameters. In the applications we assume that 
the researcher has some knowledge about at least some of the parameters 
which often represent characteristics of interest describing the process. The 
Bayesian approach incorporates this information through prior distribution 
into observed data to obtain its posterior distribution. While computation 
of the posterior distribution involves solving multidimensional integrals, the 
introduction of Markov chain Monte Carlo (MCMC) methods such as the 
Gibbs sampler and Metropolis-Hastings algorithm opened the way to analy- 
sis of complex models through decomposition and sampling from full condi- 
tional distributions; see Gamerman (1997) and Gilks et al. (1995) for general 
theory and implementation details. Some more specific discussion of the 
Bayesian dynamic modeling approach, including the choice of the hyper- 
parameters, the iterative MCMC algorithm and the implementation of the 
MCMC procedures can be found in publications by Huang and Wu (2006) 
and Wakefield (1996). The Bayesian approach was developed and tailored 
as required by the unique features of the proposed HIV dynamic models. 
The basic principles of these proposed methodologies were well established 
in the statistical literature [Gamerman (1997); Gilks et al. (1995); Wakefield 
(1996)], but the applications of these methods in this paper are nonetheless 
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innovative within the context of a system of nonhnear ODE of time- varying 
coefficient, but without a closed-form solution. 

The progress in Bayesian posterior computation due to MCMC procedures 
has made it possible to fit increasingly complex statistical models [Huang 
and Wu (2006); Wakefield (1996)] and entailed the wish to determine the 
best fitting model in a class of candidates. Thus, it has become more and 
more important to develop efficient model selection criteria. A recent publi- 
cation by Spiegelhalter et al. (2002) suggested a generalization of the Akaike 
information criterion (AIC) [Akaike (1973)] and related also to the Bayesian 
information criterion (BIC) [Schwarz (1978)] that is deviance information 
criterion (DIG). In this paper we demonstrate its usefulness to compare 
BNLME models for longitudinal HIV dynamics discussed previously. For 
completeness, a brief summary of DIG follows. More detailed discussion of 
DIG and its properties can be found in publications by Spiegelhalter et al. 
(2002) and Zhu and Garlin (2000). 

Assume that the distribution of the data, Y, depends on the parameter 
vector 'I'. Most recently, Spiegelhalter et al. (2002) suggested examining the 
posterior distribution of the deviance statistics defined by 



for Bayesian model comparison, where p(Y|^) is the likelihood function, 
that is, the conditional joint probability density function of the observed 
data Y given the parameter vector and g(y) denotes a fully specified 
standardizing term that is a function of the data alone (which thus has no 
impact on model selection). Based on the posterior distribution of D(^'), 
DIG consists of two components as follows: 



where D = £;^|y[-D(*)] = ^*|Y[-21ogp(Y|*)] and = D - L»(*) is the 
effective number of parameters, defined as the difference between the poste- 
rior mean of the deviance and the deviance evaluated at the posterior mean 
^ of the parameters. As with other model selection criteria, we caution 
that DIG is not intended for identification of the 'correct' model, but rather 
merely as a method of comparing a collection of alternative formulations. 
In our model with different baseline characteristics and/or the MEMS ad- 
herence summary metrics, DIG can be used to identify the most significant 
covariate and MEMS adherence summary metrics in contribution to viro- 
logic response. Under the model (5), in the absence of any standardizing 
function g(y), the deviance is 



21ogp(Y|*) + 21og5(Y) 



(8) 



Die = D + pD = 2D- D{^) 



2,0) = -21ogp(Y]a-2,0) 



(9) 



n 



n 



Y.<^'\yi - {^{Oi)f{yi - f.(0,)) - 



i=l 



i=l 
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As discussed above, our MCMC approach to estimating DIG first draws 
{{a~^ values from the posterior distribution, and then calculates 

corresponding {D^^^}^^^ values from (9), where G is the number of samples 
of posterior distribution. Finally, we estimate DIG as 2D — D{a~^ , 0), where 

D = h EU D^'\ = h E?=i(^-')(^\ = h Ef=i and = i = 
l,...,n}. 

3. Analysis of AIDS clinical data. 

3.1. Motivating application and observed data. The subject sample in 
our analysis was drawn from the AIDS Glinical Trials Group (AGTG) 398 
study, a randomized, double-blind, placebo-controlled, 4-Arm trial study of 
amprenavir (APV) as part of several dual protease inhibitor (PI) regimens in 
subjects with HIV infection in whom initial PI therapy had failed. Subjects 
in all arms received APV (PI), three reverse transcriptase inhibitors (RTI): 
efavirenz (EFV), abacavir (ABG) and adefovir dipivoxil (ADV) plus a second 
PI or placebo: Arm A (saquinavir = SQV), Arm B (indinavir = IDV), Arm G 
(nelfinavir = NFV) and Arm D (placebo matched for one of these three Pis). 
Subjects are HIV-infected individuals with prior exposure to approved Pis 
and who have exhibited loss of virologic suppression as reflected by a plasma 
HIV-1 RNA concentration of >1000 copies/ml. Subjects were scheduled for 
follow-up visits at study (day 0), at weeks 2, 4, 8, 12, 16 and every 8 weeks 
thereafter until week 72, and at the time of confirmed virologic failure. More 
detailed descriptions of this study and data are given by Hammer et al. 
(2002) and Pfister et al. (2003). 

As indicated previously, the primary objective of this paper is to inves- 
tigate the effect of adherence interaction with drug resistance and baseline 
covariates to prescribed ARV therapy on virologic response measured repeat- 
edly over time in HIV-infected patients. We construct a novel HIV dynamic 
model (which is parameter identifiable) with consideration of drug adher- 
ence assessed by use of MEMS data, drug susceptibility (/C50) and baseline 
covariates to link plasma drug concentration to the long-term changes in 
HIV-1 RNA observation after initiation of therapy. In the model we incor- 
porate the two clinical factors (drug adherence measured by MEMS data 
and drug susceptibility) and baseline viral load and GD4 cell count into a 
function of treatment efficacy (see Section 2.2). 

Because phenotype sensitivity testing was performed only on a subset of 
randomly selected subjects, the number of subjects available for our analysis 
was greatly reduced. We chose to consider only the subjects within Arm G 
for our analysis because this arm afforded the greatest number of subjects 
(n = 31) with available phenotypic drug susceptibility data on the two Pis 
(APV and NFV) and had available MEMS adherence data, as required for 
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our model. A summary of measurements of data to be used in our analysis 
is briefly described below. 

Plasma viral load: Plasma viral load was measured in copies/ml at de- 
signed study time by the ultrasensitive reverse transcriptase-polymerase 
chain reaction HIV-1 RNA assay (Roche Molecular Systems). Only mea- 
surements taken while on protocol-defined treatment were used in the anal- 
ysis. The exact day of viral load measurement (not predefined study week) 
was used to compute study day in our analysis. A log^o transformation was 
used in the analysis of viral load data. The graph in Figure 1 (up-left panel) 
shows the viral load trajectories of those subjects. Note that some viral load 
measurements at designed study time were not observed due to laboratory 
and other problems (for example, viral load measurement was not observed 
at week 12 for the subject displayed in Figure 1). 

Medication adherence: Medication adherence was measured by two meth- 
ods: by the use of questionnaires and by the use of MEMS [Pfister et al. 
(2003)]. Subjects completed an adherence questionnaire at study weeks. The 
questionnaire was completed by the study participant and/or by a face-to- 
face interview with study personnel. For MEMS, an MEMS cap was used to 
monitor APV and EFV compliance only. The subjects were asked to bring 
their medication bottles and caps to the clinic at each study visit, where cap 
data were downloaded to computer files and stored for later analysis. The 
MEMS adherence rate for APV was determined as the sum of positive dosing 
events divided by the sum of prescribed dosing events during the specified 
time interval. In our analysis, we assumed that NFV had the same MEMS 
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Fig. 1. The profiles of viral load measurements (in \ogiQ scale) from the 31 patients 
(up -left panel) and one trajectory of viral load (solid curve) and associated adherence rates 
(stairsteps) over time from the thirteen summary measures of MEMS data with the APV 
drug for one subject from ACTG398. 
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adherence rate as APV since both APV and NFV were prescribed with the 
same dosing schedule (twice daily), a prescribed AM and PM dosing period 
was defined for each subject and, hence, the bottles were opened twice per 
day [Pfister et al. (2003)]. As discussed previously, this study focuses mainly 
on investigating optimal strategy to summarize adherence rates determined 
by MEMS data for efficient production of virologic responses. For the MEMS 
data analysis, it was not possible to model daily adherence rates and instead 
the adherence rate was computed with the following scenarios to consider 
efi'ects of both interval length and time frame (delay of timing) for MEMS 
assessment. 

To determine the best summary metric of the MEMS adherence rate, we 
evaluated different assessment interval lengths (averaging adherence dosing 
events over 1, 2 or 3 week intervals) and different assessment time frames 
(fixing the assessment interval times to end either immediately or 1, 2 or 
3 weeks prior to the next measured viral load). Table 1 summarizes the 
MEMS assessment interval notation and definitions for the 13 scenarios. As 
an example, M2.2 in Table 1 denotes an MEMS adherence interval length 
of 2 weeks fixed to end 2 weeks prior to the next viral load measurement; 
for instance, the MEMS adherence rate for a subject at study week 8 (day 
56) was calculated as the number of nominal dosing events divided by the 
number of prescribed dosing events over study days 29-42 and this value was 
used to represent adherence from the previous study visit to the study visit 
at the day 56 for modeling. The case M serves as a reference and averages all 
the available MEMS data between viral load measurements. As an example, 



Table 1 

Summary of the MEMS interval definitions and other information 







Adherence interval definition 




Case 


MEMS adherence 


Time frame length 


Interval 


Example for week 8 




interval name 


(weeks prior to viral load 


length 


(day 56), adherence 






measurement) 




computed over days 


1 


M 


week 


visit time 


28-55 


2 


MO.l 


week 


1 week 


49-55 


3 


MO. 2 


week 


2 weeks 


42-55 


4 


MO. 3 


week 


3 weeks 


35-55 


5 


Ml.l 


1 week 


1 week 


43-49 


6 


Ml. 2 


1 week 


2 weeks 


36-49 


7 


Ml. 3 


1 week 


3 weeks 


29-49 


8 


M2.1 


2 weeks 


1 week 


36-42 


9 


M2.2 


2 weeks 


2 weeks 


29-42 


10 


M2.3 


2 weeks 


3 weeks 


22-42 


11 


M3.1 


3 weeks 


1 week 


29-35 


12 


M3.2 


3 weeks 


2 weeks 


22-35 


13 


M3.3 


3 weeks 


3 weeks 


15-35 
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Baseline Mean=29.8, SD=31.5, CV=105.8% 



Baseline Mean=61.7, SD=67.6, CV=109.5% 
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Fig. 2. The baseline (o ) and failure time (x) IC50 for APV/NFV drugs (upper panel) 
with baseline IC5Q mean, standard deviation (SD) and coefficient of variation (CV), re- 
spectively, and the baseline viral load in logj^g scale and baseline C'D4 cell count (lower 
panel) with mean, SD and CV, respectively, for the 31 subject- specific individuals from the 
ACTG398 study. Note that for a subject, if a single measurement of IC50 is observed at 
baseline only, there is no x sign appearing in the upper panel of the plot. 

the viral load (in log^g scale) and adherence rates over time from the thirteen 
cases of MEMS data with APV drug for the one representative subject are 
presented in Figure 1. 

Phenotypic virus susceptibility to drug: The phenotypic virus resistance to 
drug were retrospectively determined from baseline samples. Patients were 
selected to have samples assayed based on receiving study treatment for 
at least 8 weeks and having available sample. Some patients had virologic 
failure and phenotypic susceptibility testing done on samples at the time 
of failure. Testing was done via the recombinant virus assay (PhenoSense, 
ViroLogic Inc., South San Francisco, CA). For analysis, we used the phe- 
notype marker, /C50 [Molla et al. (1996)], to quantify agent-specific drug 
resistance. We refer to this marker as the median inhibitory concentration. 
The baseline (o) and failure time (x) /Cso's of 31 subject-specific individu- 
als for the APV and NFV drugs are displayed in Figure 2 (upper panel) and 
are used to construct IC^Q(t). Note that some subjects have only baseline 
/C50 due to the fact that they maintained viral suppression or dropped out 
from the study. If no /C50 measurement is observed at failure time for a 
subject, /C5o(t) becomes a constant in this case. 
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Baseline characteristics: The baseline viral load in log^g scale (VL) and 
the baseline CD4 cell count were chosen as covariates in the model for data 
analysis. The log-transformation of viral load is used to stabilize the variance 
of measurement error and estimation algorithm. The baseline characteristics 
of 31 subject-specific individuals with mean, standard deviation (SD) and 
coefficient of variation (CV) are displayed in Figure 2 (lower panel). To avoid 
very small (large) estimates which may be unstable, we standardized these 
covariate values. For baseline log;^Q(RNA), for instance, each logxo(RNA) 
value is subtracted by mean (4.71) and divided by standard deviation (0.70). 

3.2. Model fitting and parameter estimation results. In this section we 
apply the BNLME modeling approach to fit the data described in Section 
3.1. Based on the discussion in Section 2, the prior distribution for /x was 
assumed to be A^(J7, A) with A being a diagonal matrix. Following the idea 
of Huang and Wu (2006) for prior construction, as an example we discuss 
the prior construction for log 6. The prior constructions for other parameters 
are similar and so are omitted here. 

Ho et al. (1995) reported viral dynamic data on 20 patients; the logarithm 
of the average death rate of infected cells (log 5) is —1.125. Wei et al. (1995) 
used two different models with a group of 22 subjects to estimate death 
rate of infected cells and obtained logS with —0.84 and —1.33, respectively. 
Following these two studies, Nowak et al. (1995) estimated log 5 = —0.934 
based on 11 subjects with one possible outlying subject excluded. It can 
be seen that four estimates of log5 from these studies are —1.125, —0.84, 
—1.33 and —0.934, respectively. The individual estimates of log 5 from these 
studies approximately follow a symmetric normal distribution. Thus, we 
chose a normal distribution AA( — 1.0, 100.0) as the prior for log (5 (the large 
variance indicated that we used a noninformative prior for log 5). Similarly, 
the values of the hyper-parameters at population level are chosen as follows 
[Ho et al. (1995); Nowak et al. (1995); Nowak and May (2000); Perelson 
et al. (1996, 1997); Perelson and Nelson (1999); Verotta (2005); Wei et al. 
(1995)]: 

a = 4.5, 6 = 9.0, i/ = 10.0, 

A = diag(100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0), 

ri = (1.1, -1.0, -2.5, 1.2, 1.0, 1.0, 0.5, 0.5)^, 

n = diag(2.5, 2.5, 2.5, 2.5, 2.5, 2.5). 

We decide that one long chain is run for MCMC implication with con- 
siderations of the following two issues: (i) a number of initial "burn- in" 
simulations are discarded, since from an arbitrary starting point it would 
be unlikely that the initial simulations came from the stationary distribu- 
tion targeted by the Markov chain; (ii) one may only save every kth {k 
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being an integer) simulation sample to reduce the dependence among sam- 
ples used for parameter estimation. Because the antiviral response modeling 
involves numerically solving nonlinear differential equations, thus compu- 
tational burdens would be more pronounced with the Bayesian approaches 
via MCMC procedure. Utilizing efficient computer algorithms are critical in 
this regard. Therefore, we are going to adopt these strategies in our MCMC 
implementation using FORTRAN code that calls a differential equation sub- 
routine solver (DIVPRK) in the IMSL library (1994), which uses the Runge- 
Kutta-Verner fifth-order method. The computer codes are available from 
the corresponding author upon request. An informal check of convergence 
is conducted based on graphical techniques according to the suggestion of 
Gelfand and Smith (1990). Based on the results, we propose that, after an 
initial number of 20,000 burn-in iterations, every 4th MCMC sample was 
retained from the next 80,000 samples. Thus, we obtained 20,000 samples 
of targeted posterior distributions of the unknown parameters. 

We fitted the model to the data from 31 subjects discussed in Section 
3.1 using the proposed BNLME modeling approach. We incorporate the 
two clinical factors, drug adherence assessed by MEMS cap data and drug 
susceptibility (phenotype /C50 values), as well as baseline covariates into a 
function of drug efficacy. For model fitting, adherence rates were determined 
from MEMS data with 13 different scenarios. For model fitting and the 
purpose of comparisons, we set up a control model as the one without using 
any drug adherence, resistance and covariate information which corresponds 
to setting j{t) = 2/(exp(/3o) + 2) with IC^oit) = 1, A{t) = 1 and wi = W2 = 0; 
for this case, our model reverts to that discussed by Nowak and May (2000) 
and Perelson and Nelson (1999). The other 13 models are specified based on 
the combination of drug resistance (/C50), baseline covariate data and 13 
different adherence summary metrics listed in Table 1. 

In order to assess how adherence rates, determined from 13 different sce- 
narios, interacted with drug susceptibility and covariates to contribute to vi- 
rologic response, we fitted the models to all 13 scenarios as well as the control 
model and compared the fitting results. We found based on the DIC criterion 
(see Figure 3) that, overall, the model with adherence rate determined from 
MEMS dosing events, taken time frame length of 2 weeks prior to a viral 
load measurement with over either a 2 week assessment interval (M2.2) or a 
3 week assessment interval (M2.3), provided the best fits to the observed 
data, compared to the other 12 models for most subjects. The reference 
model with adherence rate averaged by all the available MEMS data be- 
tween viral load measurements gave a moderate fit to the observed data. 
We clearly see that all models fit the early viral load data well, but the 
control model, lacking factors for subject-specific drug adherence and sus- 
ceptibility as well as baseline covariates, failed to fit viral load rebounds 
and fluctuations, and provided the worst fitting results for the majority of 
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12 3 

Time frame length (weeks prior to RNA measurement) 

Fig. 3. Comparison of the DIG values for the models from 13 different determinants of 
MEMS adherence, interacted by drug resistance and covariates, with the control model. 
The two horizontal lines represent the DIG values for the control model and the reference 
model with adherence rate determined by case M, respectively. 

subjects. For the purpose of illustration, the model fitting curves from the 
control model (solid curves), the best fit model (M2.2: dotted curves) and 
the reference model (M: dashed curves) are displayed in Figure 4 for the 
three representative subjects. 

For the purpose of comparison, Figure 5 presented the population posterior 
means and the corresponding 95% equal-tail credible intervals (CI) of the eight 
parameters for the control model, the best fit model (M2.2) and the reference 
model (M). For the six dynamic parameters {c,d,dT, p, R, /3o), it is shown 
that the population estimates for the control model have higher clearance 
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Fig. 4. The estimate of viral load trajectory from the model fitting with the 3 different de- 
terminants of adherence: (i) Gontrol model (solid curves), (ii) M2.2 model (dotted curves) 
and (iii) reference (M) model (dashed curves) for the three representative subjects. The 
observed values are indicated by circles. 
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1=Control 
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Fig. 5. A summary of the estimated posterior means (o) of population parameters and 
the corresponding 95% equal-tail credible intervals ( CI) for the models from 3 different 
determinants of adherence. 

rate of free virions (c), lower death rate of infected cells {5), higher death rate 
of target T cells (dr), smaller viral load scaling factor (p), higher basic repro- 
ductive ratio for the virus {R) and larger (f) than those for the best fit model 
(M2.2) and reference model (M), while the population estimates for the best 
fit model and reference model are generally similar. These differences may 
result from the effects of drug adherence interacted with drug resistance 
and covariates in the models. For the other two covariate effect parame- 
ters (/3i,/32) which are relevant to treatment effect, we will discuss them 
separately in Section 3.4. In terms of the individual parameter estimates, a 
large between-subject variation in the estimates of all individual dynamic 
parameters was observed (data not shown here). Overall, the coefficient of 
variation ranges from 15.4% to 88.9% for all parameters. 



3.3. Effects of adherence rates determined by different MEMS summary 
metrics. Figure 3 in Section 3.2 displayed a comparison of the DIG val- 
ues for the models from 13 different determinants of MEMS adherence, in- 
teracted by drug resistance and covariates, with the control model. The 
observed patterns shown in Figure 3 provided information to answer the 
following questions: (i) what MEMS assessment interval length is best and 
(ii) what MEMS assessment time frame (delay effect of timing) is best? 
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We can see that when the time frame for MEMS assessment is fixed, mod- 
els with a 2 week MEMS assessment interval length generally outperform 
models with an assessment interval length of 1 or 3 weeks except for the 
time frame length with 3 weeks prior to viral load measurement where the 
model with a 3 week MEMS assessment interval length performs best. 

Regardless of the assessment interval length, models which assess compli- 
ance 2 weeks prior to viral load measurement generally outperform models 
which assess compliance immediately before viral load measurement, 1 week 
before or 3 weeks before viral load measurement. Overall, the model with a 
MEMS assessment interval length of 2 weeks measured from 4 to 2 weeks 
prior to viral load measurement (M2.2) was significantly a better predicator 
of viral load over time than any other models, with the exception of the 
M2.3 model which shows no significant difference from the M2.2 model in 
terms of DIG values. 

3.4. Treatment effects of baseline characteristics interacted with clinical 
factors. Figure 5 summarized the population posterior means and the cor- 
responding 95% equal-tail CI of the covariate effect parameters /3i and (^2 
for the best fit model (M2.2) and the reference model (M). It can be seen 
that estimates of /3i (coefficient of baseline viral load) are negative, while 
estimates of P2 (coefficient of baseline CD4 cell count) are positive. In fact, 
other models also provided the same scenarios for the estimates of these two 
parameters (not shown here). As an example, we report results based on the 
best fit model (M2.2). We can observe from Figure 5 that the estimates of /3i 
and P2 are /3i = -0.67 with 95% CI (-1.056, -0.193) and /32 = 0.719 with 
95% CI (0.371, 1.058). It indicates that, according to antiviral drug efficacy 
model (4), the baseline viral load (/3i = —0.67) has a significant positive 
effect on drug efficacy 7(t), while the baseline CD4 cell count (/32 = 0.719) 
has a significant negative effect on 7(t) since the corresponding 95% credible 
intervals do not contain zero for both parameters. These findings could sug- 
gest us with the following different ways. The lowest value of 7(t) [highest 
(j) as displayed in Figure 6(b)] occurs in the subjects with the best prognosis 
(higher baseline CD4 cell count and lower baseline viral load). Alternatively, 
the highest value of 7(t) (lowest (p) occurs in those with the worst prognosis 
(lower baseline CD4 cell count and higher baseline viral load). A possible 
explanation is that there is a fioor effect of viral load (or ceiling/fioor effect 
of CD4 cell count) that is not captured in the model. Further, given that 
baseline CD4 cell count and viral load are jointly used to make treatment 
decisions and are known to be negatively correlated as shown in Figure 6(a), 
the result based on the combination of baseline viral load and CD4 count 
in 7(t) indicates that the baseline CD4 cell count and viral load have the 
opposite effect on drug efficacy which might be intuitively understandable. 
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Fig. 6. (a) Correlation between (standardized) baseline logjQ (RNA) and CD4 cell count. 
The correlation coefficient (r) and p-value are obtained from the Spearman rank correlation 
test. The line is a robust (MM-estimator) linear regression fit. (b) Estimated (j) values for 
the 31 subject- specific individuals. 

4. A simulation study. In this paper we investigated the association be- 
tween virologic outcomes and medication adherence with confounding fac- 
tors based on the data from 31 subjects. As both one referee and the asso- 
ciate editor suggested that a simulation study may be useful to evaluate how 
our method performs, in this section we conduct a limited simulation study 
here due to intensive computations involved. The scenario we consider is as 
follows. 

We simulate a clinical trial with 31 patients receiving antiviral treatment. 
For each patient, we assume that the designs of this experiment, in particu- 
lar, the sampling times for viral load, were the same as those in the ACTG 
398 study. The data for the phenotype marker (baseline and failure /Cso's), 
medication adherence and the baseline viral load/CD4 cell count were taken 
from the ACTG 398 study, where medication adherence was calculated by 
the M2.2 summary measure. The "true" values of unknown parameters were 
the same as those estimated from the data set of 31 subjects which were re- 
ported in Section 3. With generated individual true parameters based on the 
equation (6), known data [IC^o{t),A{t).,wi and W2\-, we generated random 
samples for response (viral load) based on model (2). The values of hyper- 
parameters are chosen to be the same as those in Section 3. For each simu- 
lated data set, we fit the model using the Bayesian approach. The MCMC 
techniques consisting of a series of Gibbs sampling and M-H algorithms were 
the same as those in the real data analysis. We performed 50 replications 
and obtained the mean estimates (ME) of population parameters together 
with the corresponding relative bias (RB), which is the difference between 
the mean estimate and the true value of the parameter divided by absolute 
value of the true parameter, and the standard error (SE), defined as the 
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The true values (TV) of parameters and mean estimates (ME) of population parameters 
with 50 replications as well as the corresponding relative bias (RB), defined as 
100 X [ME - TV)/\ TV\, and standard error (SE), defined as 100 x VMSE/\ TV\ 



Parameter 


TV 


ME 


RB (%) 


SE (%) 


logc 


0.767 


0.771 


0.522 


9.127 


log 5 


-0.977 


-1.013 


-3.685 


7.971 


log dr 


-4.086 


-4.101 


-0.367 


5.871 


logp 


0.433 


0.388 


-10.39 


18.03 


logi? 


1.040 


1.100 


5.769 


3.089 




-2.615 


-2.604 


0.421 


4.902 


Pi 


-0.670 


-0.665 


0.746 


10.13 


P2 


0.719 


0.698 


-2.921 


13.09 



square root of mean-squared error divided by the absolute value of the true 
parameter. 

In Table 2 we summarize the true values (TV) of parameters and the ME 
of population parameters with 50 replications as well as the correspond- 
ing RB and SE. The percentage is based on the absolute value of the true 
parameter. It can be seen from Table 2 that the RB (%) for population 
parameter estimates are very small, ranging from 0.522 to 10.39, and the 
SE (%) ranges from 3.089 to 18.03. The simulation results indicate that our 
method with considering the M2.2 model performs reasonably well in terms 
of estimates of parameters except for the viral load proportionality factor 
log/j which has larger RB and SE. That is, our method produces a substan- 
tially biased estimate and may severely underestimate log p. This may be 
explained by the fact that it is probably caused from inaccurate numerical 
solutions to the system of ODE (2) which was used to construct the BNLME 
model. 

5. Concluding discussion. In developing long-term dynamic modeling, 
this paper introduced a dynamic mechanism specified by a system of time- 
varying ODE to (i) establish a link between success of ARV therapy in 
virologic response and MEMS adherence confounded by drug resistance and 
baseline covariates, (ii) fully integrate viral load, MEMS adherence, drug re- 
sistance and baseline covariates data into the statistical inference and analy- 
sis, and (iii) provide a powerful tool to evaluate the effects of MEMS adher- 
ence determined by a different summary metric on virologic response using 
the BNLME modeling approach. This approach cannot only combine prior 
information with current clinical data for estimating dynamic parameters, 
but also deal with complex dynamic systems. Thus, the results of estimated 
dynamic parameters based on this model should be more reliable and rea- 
sonable to interpret long-term HIV dynamics. Our models are simplified 
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with the main goals of retaining crucial features of HIV dynamics and, at 
the same time, guaranteeing their applicability to typical clinical data, in 
particular, long-term viral load measurements. The proposed model fitted 
the clinical data reasonably well for most patients in our study, although the 
fitting for a few patients was not completely satisfactory because of unusual 
viral load fiuctuation patterns for these subjects. 

We have explored the practical performance of DIG for the comparison 
of developed models. DIG, a Bayesian version of the classical deviance for 
model assessment, is particularly suited to compare Bayesian models whose 
posterior distribution has been obtained using MGMG procedures and can 
be used in complex hierarchical models where the number of unknowns often 
exceeds the number of observations and the number of free parameters is not 
well defined. This is in contrast to AIG and BIG, where the number of free 
parameters needs to be specified [Zhu and Garlin (2000)]. Overall, combined 
with more traditional residual analysis and posterior predictive model checks 
as discussed in this paper, DIG appears to offer a comprehensive framework 
for comparison and evaluation within a complex model class. 

Several studies investigated the association between virologic responses 
and adherence assessed by MEMS data only without considering other con- 
founding factors such as drug resistance using standard modeling methods 
including Poisson regression [Knafi et al. (2004)], logistic regression [Vrijens 
et al. (2005)] and the linear mixed-effects model [Liu et al. (2007)]. In this 
article we employed the proposed dynamic model and associated BNLME 
modeling approach to assessment of effects of adherence determinants based 
on MEMS dosing events in predicting virologic response. In particular, we 
investigated (i) how to summarize the MEMS adherence data for efficient 
prediction of virological response after accounting for potential confound- 
ing factors such as drug resistance and baseline covariates, and (ii) how to 
evaluate treatment effect of baseline characteristics interacted with MEMS 
adherence and other clinical factors. Note that a further study in comparing 
the performance of these different methods may be important and war- 
ranted, although some challenges are observed in terms of different model 
structures and data characteristics. 

The results indicate that the best summary metric for prediction of viro- 
logic response based on DIG criterion is the adherence rate determined by 
MEMS dosing events averaged over an assessment interval of 2 or 3 weeks, 
and 2 weeks prior to the next measured viral load observation (denoted by 
M2.2 or M2.3). We found that the best MEMS adherence predictor (M2.2) 
of the effectiveness of ARV medications on virologic response is consistent 
with that reported in Huang et al. (2008) in which, however, the next best 
MEMS adherence predictor (Ml. 2) is different from what is obtained in this 
paper. This difference may be due to the various reasons as follows. In the 
study by Huang et al. (2008), (i) it directly applied the model (1) to fit data 
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and, thus, some assumptions were made due to parameter unidentifiable is- 
sues; (ii) the analysis used the mean of the sum of squared deviations as a 
criterion to evaluate model fitting results; (iii) it assumed /C50 data were 
extrapolated linearly to the whole treatment period instead of a log-linear 
extrapolation offered in this paper which is considered more reasonable bi- 
ologically; and (iv) it did not incorporate baseline covariates in the model. 
In addition, the superiority of the M2.2 model, associated with the MEMS 
adherence rate based on time frame length of 2 weeks prior to a viral load 
measurement with over a 2 week assessment interval, may be explained by 
the fact that it probably reflects how long it takes for resistance mutations 
to first arise and then come to dominate the plasma population of a virus. As 
pointed out by an anonymous referee, this finding may also be interpreted 
as follows. Low adherence two weeks prior to the viral load measurement 
may not have had sufficient time for viral rebound to occur. 

In this paper we set up a connection between subject-specific baseline 
characteristics with interaction of clinical factors and drug efficacy. We also 
found that, according to antiviral drug efficacy model (4), the baseline viral 
load had a positive effect on drug efficacy, while the baseline CD4 cell count 
had a negative effect on it. Our results may be explained by the fact that 
for those patients with higher baseline viral load, the drug efficacy needs 
to be higher than that for those with lower baseline viral load. Therefore, 
a strong treatment is recommended for those patients with higher baseline 
viral load. On the other hand, patients with higher CD4 cell count may need 
lower drug efficacy so that a more potent ARV drug regimen is not necessary 
for these patients to avoid side-effects of drug use. The results may suggest 
the benefit of initiating ARV therapy with a lower baseline viral load and/or 
a higher baseline CD4 count. These results coincide with those investigated 
by Notermans et al. (1998) and Wu et al. (2005) whose results were ob- 
tained using correlation analysis. Note that given the estimated parameters, 
the subject with both a high baseline CD4 cell count and a relatively high 
baseline viral load [upper right quadrant of Figure 6(a)] has a very different 
(j) than that with a similar baseline CD4 cell count, but a low baseline viral 
load [upper left quadrant of Figure 6(a)]. It is possible that the subject in 
the upper right quadrant was more recently infected (hence the higher base- 
line CD4 cell count) or perhaps with a drug resistant virus and would not 
be a candidate for a regimen with a "less potent drug efficacy." 

Our findings need to be interpreted in light of the study limitations. First, 
in the ACTG 398 study, because phenotype sensitivity testing was performed 
only on a subset of randomly selected subjects, we chose 31 patients who 
have available data for analysis in this paper. Second, due to reasons such as 
lost caps and malfunction of caps, there were inaccurate MEMS data across 
the treatment period which may not reflect actual adherence profile for in- 
dividual patients and, thus, the data quality could have some impact on the 
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results. Third, because of technical limitations, the undetectable values of vi- 
ral load were replaced with 25 copies/ml for analyses, which could introduce 
some bias due to a cluster of ties of data points. Finally, this paper combined 
new technologies in mathematical modeling and statistical inference with ad- 
vances in HIV/ AIDS dynamics and ARV therapy to quantify complex HIV 
disease mechanisms. The complex nature of HIV/ AIDS ARV therapy will 
naturally pose some challenges including missing data and measurement er- 
ror in clinical factors and covariates. These complicated problems, which are 
beyond the scope of this article, may be addressed, for example, using the 
joint model method [Wu (2002)] and other techniques [Carroll et al. (1995)], 
and are warranted for further investigation. Nevertheless, these limitations 
would not offset the major findings from this study. 

As the Associate Editor pointed out, we assumed that the distributions 
of the random error and random effects are normal, which is a common as- 
sumption in the literature for statistical inference. However, due to the na- 
ture of AIDS clinical data, it is possible that the data may contain outliers 
and/or depart from normality and, thus, statistical inference and analysis 
with normal assumption may lead to misleading results [Verbeke and Lesaf- 
fre (1996); Ghosh et al. (2007)]. Specially non-normal characteristics such as 
skewness with heavy right or left tail may appear often in virologic responses. 
Thus, a normality assumption may be too restrictive to provide an accurate 
representation of the structure that is often present in repeated measures 
and clustered data. Thus, it is of practical interest to investigate nonlinear 
models with a skew-normal distribution or t distribution for (within-subject) 
random error and random effects which are more robust to outliers and skew- 
ness than those with a normal distribution. In our recent study [Huang and 
Dagne (2010)] we addressed a Bayesian approach to nonlinear mixed-efi'ects 
models in conjunction with the HIV dynamic model and relaxed the nor- 
mality assumption by considering both random error and random-effects to 
have a multivariate skew-normal distribution. The proposed model provides 
flexibility in capturing a broad range of non-normal behavior and includes 
normality as a special case. The results suggest that it is very important to 
assume models with a skew-normal distribution in order to achieve robust 
and reliable results, in particular, when the data exhibit skewness. We are 
actively applying this methodology into the data investigated in this paper 
and will report the results in a future study. 

In summary, the mechanism-based dynamic model is powerful and effi- 
cient to characterize relations between antiviral response and medication 
adherence, drug susceptibility as well as baseline characteristics, although 
some biological assumptions are required. It is important to find a way to 
incorporate subject-specific information with regard to drug susceptibility, 
medication adherence and baseline characteristics in predicting long-term 
virologic response. Since each of these factors may only contribute a very 



DYNAMIC BAYESIAN NONLINEAR MIXED-EFFECTS MODEL 



25 



small portion to virologic response and they may be confounded through 
complicated interactions, the appropriate modeling of the combination ef- 
fects of these factors is critical to efficiently utilize the information in viro- 
logic response predictions. The viral dynamic model and associated statis- 
tical approaches discussed here provide a good avenue to fulfill this goal. In 
particular, MEMS adherence rate summarized by an optimal way in terms 
of assessing both interval lengths and time frame lengths prior to viral load 
measurement is an important factor that significantly determines the effec- 
tiveness of ARV treatment and needs to be taken into account in analy- 
sis of virologic responses. Our results demonstrate that MEMS adherence 
data may not predict virologic response well unless the MEMS cap data are 
summarized in an appropriate way as reported in Section 3.3. Additionally, 
although this paper concentrated on HIV dynamics, the basic concept of 
longitudinal dynamic systems and the proposed methodologies in this paper 
are generally applicable to dynamic systems in other fields such as biology, 
medicine, engineering or PK/PD studies as long as they meet the relevant 
technical specification — a system of ODE. 
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